Generalized Linear Models - Logistic Regression¶
Odds, Log Odds, Logit, Logistic¶
Regression¶
Link¶
- We need to "link" our regression to the bernoulli family
- We want to predict a 0 or 1 for each class
- Y is the probability of being in a particular class
- Our original regression doesn't satisfy
- Binary data doesn't follow the normal distribution
- Probabilities usually aren't linear
- Our values extends beyond 0 and 1
Odds¶
Logistic¶
Sigmoid¶
Logistic Regression¶
Setup¶
In [118]:
import pandas as pd
import numpy as np
import pandas_profiling
import matplotlib.pyplot as plt
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.feature_selection import RFE
from sklearn.metrics import auc
from sklearn.metrics import f1_score
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from sklearn.metrics import roc_curve
from sklearn.metrics import roc_auc_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import log_loss
from sklearn.metrics import cohen_kappa_score
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import OneHotEncoder
from sklearn.impute import SimpleImputer
import statsmodels.api as sm
from yellowbrick.classifier import ROCAUC
from yellowbrick.classifier import ConfusionMatrix
from yellowbrick.classifier import ClassificationReport
from yellowbrick.classifier import ClassPredictionError
from yellowbrick.classifier import DiscriminationThreshold
from yellowbrick.features import Manifold
Data¶
In [48]:
breastCancer = load_breast_cancer()
print(breastCancer.DESCR)
In [49]:
breastCancerDF = pd.DataFrame(breastCancer['data'], columns=breastCancer['feature_names'])
breastCancerDF['TARGET'] = breastCancer['target']
breastCancerDF.head()
Out[49]:
EDA¶
In [50]:
pandas_profiling.ProfileReport(breastCancerDF)
Out[50]:
Visualization¶
T-SNE¶
- https://mlexplained.com/2018/09/14/paper-dissected-visualizing-data-using-t-sne-explained/
- t-SNE is not a linear projection (aka PCA)
- t-SNE creates a low-dimensional mapping allowing it to capture non-linearity
- Used here to get some intuition about the data large dimension space
In [51]:
X = breastCancerDF.drop('TARGET', axis=1)
y = breastCancerDF['TARGET']
plt.figure(figsize=(15,10))
visualizer = Manifold(manifold='tsne', target='discrete')
visualizer.fit_transform(X, y)
visualizer.poof()
Classification¶
Correlated Features¶
- Remove highly correlated features identified in pandas_profiling
In [52]:
select = [
#'mean_radius',
'mean_texture',
#'mean_perimeter',
'mean_area',
'mean_smoothness',
'mean_compactness',
#'mean_concavity',
'mean_concave_points',
'mean_symmetry',
'mean_fractal_dimension',
#'radius_error',
'texture_error',
#'perimeter_error',
'area_error',
'smoothness_error',
'compactness_error',
'concavity_error',
'concave_points_error',
'symmetry_error',
'fractal_dimension_error',
#'worst_radius',
#'worst_texture',
#'worst_perimeter',
#'worst_area',
'worst_smoothness',
'worst_compactness',
'worst_concavity',
#'worst_concave_points',
'worst_symmetry',
'worst_fractal_dimension'
]
Test and Train Split¶
- Cut into Train & Test -- stratify Target
In [53]:
X = breastCancerDF[select]
y = breastCancerDF['TARGET']
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=1337)
Pipeline¶
- No need to scale variables for Logistic Regression
In [54]:
featuresNumeric = X.columns.tolist()
transformerNumeric = Pipeline(steps=[
('imputer', SimpleImputer(strategy='median'))
])
featuresCategorical = []
transformerCategorical = Pipeline(steps=[
('imputer', SimpleImputer(strategy='constant', fill_value='missing')),
('onehot', OneHotEncoder(handle_unknown='ignore'))
])
preprocessor = ColumnTransformer(transformers=[
('numeric', transformerNumeric, featuresNumeric),
('categorical', transformerCategorical, featuresCategorical)
])
GLM¶
Model¶
In [55]:
logistic = LogisticRegression()
logPipe = Pipeline([
('preprocess', preprocessor),
('logistic', logistic)
])
logPipe.fit(X_train, y_train)
Out[55]:
Predictions¶
In [56]:
preds = logPipe.predict(X_test)
proba = logPipe.predict_proba(X_test)
print('Score Train {:.3f}'.format(logPipe.score(X_train, y_train)))
print('Score Test {:.3f}'.format(logPipe.score(X_test, y_test)))
In [57]:
predDF = pd.DataFrame(proba, columns=['Probability False', 'Probability True'])
predDF['Predicted Class'] = preds
predDF.head()
Out[57]:
Threshold¶
- Defaulted to 0.5, but can be set from 0.0 to 1.0 depending on "cost" of error
- Precision: An increase in precision is a reduction in the number of false positives
- Recall: An increase in recall decreases the likelihood that we miss a positive class
- F1 Score: The F1 score is the harmonic mean between precision and recall
- Queue Rate: This metric describes the percentage of instances that must be reviewed.
In [58]:
plt.figure(figsize=(15,10))
visualizer = DiscriminationThreshold(logPipe)
visualizer.fit(X_train, y_train) # Fit the training data to the visualizer
visualizer.poof() # Draw/show/poof the data
Errors¶
Accuracy¶
- Proportion of labels we correctly identify
In [92]:
print('Test Accuracy = {:.3f}'.format(logPipe.score(X_test, y_test)))
Log Loss¶
- Common loss objective function
- Quantifies the accuracy of a classifier by penalising false classifications
- Minimising the Log Loss == Maximising the classifier accuracy
- Log Loss penalizes erroneous predictions made with high confidence more than those with low confidence
- Class = 1 & Predict Prob Class 1 = 0.5, Log Loss == 1
- Class = 1 & Predict Prob Class 1 = 0.0, Log Loss == Very Large (+Infinity)
- Class = 1 & Predict Prob Class 1 = 1.0, Log Loss == Very Small (Approaching 0)
In [59]:
print('Log Loss: {:.3f}'.format(log_loss(y_test, preds)))
Confusion Matrix¶
In [60]:
confusion_matrix = confusion_matrix(y_test, preds)
print(confusion_matrix)
Classification Report¶
In [61]:
rpt = classification_report(y_test, preds)
print(rpt)
AUC - ROC Curve¶
- https://towardsdatascience.com/understanding-auc-roc-curve-68b2303cc9c5
- ROC curve measures performance on classification problem at various thresholds
- ROC is a probability curve plotted with TPR against FPR
- AUC represents degree or measure of separability
- Total Area = 1.0
- Random = 0.5
- Higher AUC equates to a better the model
In [62]:
logit_roc_auc = roc_auc_score(y_test, preds)
fpr, tpr, thresholds = roc_curve(y_test, proba[:,1])
plt.figure(figsize=(15,10))
plt.plot(fpr, tpr, label='Logistic Regression (area = %0.2f)' % logit_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([-0.01, 1.0])
plt.ylim([0.0, 1.01])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC Curve')
plt.legend(loc="lower right")
plt.style.use('seaborn-whitegrid')
plt.show()
AUC¶
In [63]:
print('AUC: {:.3f}'.format(logit_roc_auc))
Cohen's Kappa¶
- Baseline is found by always guessing the majority class
- In our example we have 90 Positives out of 143 in test (~ 63%)
- Intuition : If our accuracy were .63, our Cohen's Kappa would = 0
- We can include weights on the errors we make
In [123]:
print('Accuracy on Test : {:.3f}'.format(acc))
print('Baseline on Test : {:.3f}'.format(base))
print("Cohen's Kappa : {:.3f}".format(cohen_kappa_score(logPipe.predict(X_test), y_test)))
Yellowbrick Visualization¶
Confusion Matrix¶
In [64]:
plt.figure(figsize=(8,6))
visualizer = ConfusionMatrix(logPipe, classes=[0,1], percent=True, fontsize=20)
visualizer.fit(X_train, y_train)
visualizer.score(X_test, y_test)
visualizer.poof()
Classification Report¶
In [65]:
plt.figure(figsize=(10,8))
visualizer = ClassificationReport(logPipe, classes=[0,1], support=True)
visualizer.fit(X_train, y_train) # Fit the visualizer and the model
visualizer.score(X_test, y_test) # Evaluate the model on the test data
g = visualizer.poof() # Draw/show/poof the data
ROC Curve¶
In [66]:
plt.figure(figsize=(15,10))
visualizer = ROCAUC(logPipe, classes=['No Cancer', 'Cancer'])
visualizer.fit(X_train, y_train)
visualizer.score(X_test, y_test)
g = visualizer.poof()
Class Prediction Error¶
In [67]:
plt.figure(figsize=(8,6))
visualizer = ClassPredictionError(logPipe, classes=[0,1])
visualizer.fit(X_train, y_train)
visualizer.score(X_test, y_test)
g = visualizer.poof()
Interpretation¶
20 Variable Model¶
In [68]:
model = sm.Logit(y_train, X_train)
result = model.fit()
print(result.summary2())
preds20 = result.predict(X_test)
Coefficients¶
- Coefficients affect Log Odds - Exponentiate to see change
- 1 Unit change of coefficient relates to Log-Odds of Cancer
- Coefficient changes are not linear!
- Exponentiate coefficient to find the change in Odds
- Slope is instructive (Negaitve is a reduction in odds)
Goodness of Fit¶
- Deviance
- Hosmer–Lemeshow test
- Wald statistic
- Likelihood ratio test
Observation¶
- There are several extremely large coefficients and serveral with large p-values (intervals including 0).
- fractal_dimension_error 2,434
- smoothness_error 819
- mean_fractal_dimension 673
- We could be overfitting or including too many features in the model
- We should look at:
- Feature selection
- Regularization
Feature Selection¶
Scoring Top N Features¶
In [87]:
metricAUC = []
metricF1 = []
metricLogLoss = []
for features in range(20):
rfe = RFE(logistic, features + 1, 1)
rfe = rfe.fit(X_train, y_train)
features_bool = np.array(rfe.support_)
features = np.array(X.columns)
topFeatures = features[features_bool].tolist()
model = sm.Logit(y_train, X_train[topFeatures])
result = model.fit()
preds = result.predict(X_test[topFeatures])
fpr, tpr, thresholds = roc_curve(y_test, preds, pos_label=1)
metricAUC.append(auc(fpr, tpr))
metricF1.append(f1_score(y_test, preds > 0.5))
metricLogLoss.append(log_loss(y_test, preds))
Top N AUC¶
In [70]:
pd.DataFrame(metricAUC, columns=['Features']).plot(figsize=(15,10))
plt.title('F1 Score for Top N Features')
plt.ylabel('AUC')
plt
Out[70]:
Top N F1-Score¶
In [71]:
pd.DataFrame(metricF1, columns=['Features']).plot(figsize=(15,10))
plt.title('F1 Score for Top N Features')
plt.ylabel('F1 Score')
plt
Out[71]:
Top N Log Loss¶
In [89]:
pd.DataFrame(metricLogLoss, columns=['Features']).plot(figsize=(15,10))
plt.title('LogLoss for Top N Features')
plt.ylabel('LogLoss Score')
plt
Out[89]:
Top 15 Features¶
- Best AUC
- Lowest LogLoss
In [72]:
logistic = LogisticRegression()
logPipeRFE = Pipeline([
('preprocess', preprocessor),
('logistic', logistic)
])
In [73]:
rfe = RFE(logistic, 15, 1)
rfe = rfe.fit(X_train, y_train)
features_bool = np.array(rfe.support_)
features = np.array(X.columns)
topFeatures = features[features_bool].tolist()
topFeatures
Out[73]:
15 Variable Model¶
In [74]:
model = sm.Logit(y_train, X_train[topFeatures])
result = model.fit()
print(result.summary2())
preds15 = result.predict(X_test[topFeatures])
Comparison¶
- The 15 feature model includes several features that have a high p-value
- High p-value or intervals stretching across zero potentially could be removed
- Alternatively, we could look to regularize the coefficients
In [75]:
resultDF = pd.DataFrame(columns=['LogLoss', 'AUC'])
20 Features¶
In [76]:
fpr, tpr, thresholds = roc_curve(y_test, preds20, pos_label=1)
area_under_curve = auc(fpr, tpr)
logloss = log_loss(y_test, preds20)
resultDF.loc['20 Feature Model'] = [logloss, area_under_curve]
print('AUC: {:.3f}'.format(area_under_curve))
print('Log Loss: {:.3f}'.format(logloss))
In [77]:
preds20_class = preds20 > 0.5
print(classification_report(y_test, preds20_class))
15 Features¶
In [78]:
fpr, tpr, thresholds = roc_curve(y_test, preds15, pos_label=1)
area_under_curve = auc(fpr, tpr)
logloss = log_loss(y_test, preds15)
resultDF.loc['15 Feature Model'] = [logloss, area_under_curve]
print('AUC: {:.3f}'.format(area_under_curve))
print('Log Loss: {:.3f}'.format(logloss))
In [79]:
preds15_class = preds15 > 0.5
print(classification_report(y_test, preds15_class))
Regularization¶
Grid Search Regularization¶
In [80]:
logistic = LogisticRegression()
regPipe = Pipeline([
('preprocess', preprocessor),
('logistic', logistic)
])
param_grid = {
'logistic__penalty' : ['l1'],
'logistic__C' : np.logspace(-4, 4, 20)
}
clf = GridSearchCV(regPipe, param_grid = param_grid, cv = 5, verbose=True, n_jobs=-1)
best_clf = clf.fit(X_train, y_train)
In [81]:
best_clf.best_params_
Out[81]:
Regularized Model¶
In [82]:
model = sm.Logit(y_train, X_train)
result = model.fit_regularized(method='l1', alpha=1.0, L1_wt=29.763)
print(result.summary2())
predsReg = result.predict(X_test)
Errors¶
In [83]:
fpr, tpr, thresholds = roc_curve(y_test, predsReg, pos_label=1)
area_under_curve = auc(fpr, tpr)
logloss = log_loss(y_test, predsReg)
resultDF.loc['Regularized Model'] = [logloss, area_under_curve]
print('AUC: {:.3f}'.format(area_under_curve))
print('Log Loss: {:.3f}'.format(logloss))
In [84]:
predsReg_class = predsReg > 0.5
print(classification_report(y_test, predsReg_class))
Overall Results¶
In [86]:
resultDF.sort_values('LogLoss')
Out[86]:
Advice & Tips¶
- Advantages:
- Computationally efficient
- Interpretable
- Does not require features scaling
- Does not require tuning
- Easy to regularize
- Logistic Regression works best when:
- Remove features that are highly correlated with eachother
- Remove features that have little relationship with the dependent variable
- Downsides
- Can be easily outperformed by complex models
- Can’t solve non-linear problems since it has a linear decision surface / boundary
- Practical Tips
- Consider feature selection
- Test Regularization L1 & L2
- Use statsmodel for interpreation
In [ ]: